DIPnet Population Structure Notebook

Setup

library(seqinr)
library(ape)  
library(pegas)
library(hierfstat)
library(mmod)
library(adegenet)
library(plyr)
library(strataG)
library(iNEXT)
library(gdm)
library(gdistance)
library(ecodist)
library(dplyr)
library(reshape2)
library(WriteXLS)
library(ggplot2)
library(knitr)
library(vegan)
source("config.R")
source("DIPnet_Stats_Functions.R")
rescale<-function(x){
  normalized<-(x-min(x))/(max(x)-min(x))
  return(normalized)
}

######################################################################
# 1. Import the IPDB and Fst tables
ipdb<-read.table(ipdb_path,sep="\t",header=T,stringsAsFactors = F,quote="", na.strings=c("NA"," ","")) 
#read in geographical regionalizations from Treml
spatial<-read.table(spatial_path, header=T, sep="\t",stringsAsFactors = F, na.strings=c("NA"," ",""), quote="")
#read in geographical regionalizations from Beger
spatial2<-read.table(spatial2_path, header=T,sep="\t", stringsAsFactors = F, na.strings=c("NA"," ",""), quote="")
#read in ABGD groups
abgd<-read.table(abgd_path, header=T, sep="\t", stringsAsFactors = F)
#join spatial
ipdb<-join(ipdb,spatial, by = "IPDB_ID",type = "left")
ipdb<-join(ipdb,spatial2[,c(2,18:24)], by = "IPDB_ID", type = "left")
#join ABGD
ipdb<-join(ipdb,abgd[,c(1,3)], by = "IPDB_ID",type = "left")
# drop hybrids and divergent individuals
ipdb<-ipdb[ipdb$IPDB_ID %in% drops == FALSE, ] 
## remove anything not included in the ecoregions scheme (some dolphins, some COTS from Kingman and Madagascar(?), some A. nigros from Kiribati, som C. auriga from Fakareva, hammerheads from Western Australia, and West Africa, and some dolphins from the middle of the eastern tropical pacific
ipdb_ecoregions<-ipdb[-which(is.na(ipdb$ECOREGION)),]
## remove anything that doesn't occur in the 3 Indo-Pacific realms
ipdb_ip<-ipdb_ecoregions[which(ipdb_ecoregions$REALM %in% c("Central Indo-Pacific","Western Indo-Pacific","Eastern Indo-Pacific")),]
amova_ts_path<-"/Users/eric/google_drive/DIPnet_Gait_Lig_Bird/DIPnet_WG4_first_papers/statistics/By_Species/Hierarchical_structure"
amova_abgd_path<-"/Users/eric/google_drive/DIPnet_Gait_Lig_Bird/DIPnet_WG4_first_papers/statistics/By_ESU/Hierarchical_structure"
# Read in Veron barriers
barriers<-read.csv("VeronBarriers.csv",header=F,stringsAsFactors = F)

Introduction

This is an R Markdown Notebook. I am trying it out as a way to reproducibly document my work on the DIPnet Population Structure Paper.

The overall workflow for this analysis is as follows:

• Run AMOVA according to several different regionalizations
  • Briggs and Bowen 2013 - Fish Biogeography
  • Veron et al. 2015 - Coral Biogeography
  • Spalding et al. 2007 -
    • Marine Ecoregions
    • Provinces
    • Realms
  • Kulbicki et al. 2013 -
    • Regionalization based on all species,
    • or just “reliable” species
  • Keith et al. 2013 - Faunal Breaks
• Select the best model for each species according to the Bayesian Information Criterion
• Analyze individual “barriers” as implied by the regionalization that best explains the distribution of variation in the most species by modeling pairwise genetic distances given geographic distance and each barrier (individually or together):
  • Generalized Dissimilarity Modeling (GDM)
  • Multiple Regression with Distance Matrices (MRDM)
  • Distance Based Redundancy Analysis (dbRDA)
• Special considerations:
  • What is a species? Do the analysis using:
    • Traditional species boundaries
    • Automated Barcode Gap Discovery (ABGD) to account for cryptic species
  • Include sequence distance in AMOVA?
    • No: traditional Fst, specifically Weir and Cockerhams Theta
    • Yes: Phi-ST

Run AMOVA loops

Given the special considerations, we need to run 4 different flavors of AMOVA: 1. WCTheta and TradSpec 2. WCTheta and ABGD 3. PhiST and TradSpec 4. PhiST and ABGD

Thus, we need to loop through all of the regionalizations above, running all 4 flavors of AMOVA.

## Loop through hypotheses, calculating AMOVA
hypotheses<-c("ECOREGION", "PROVINCE","REALM","Bowen","Keith","Kulbicki_r","Kulbicki_b", "VeronDivis")
amova_list<-list()

for(h in hypotheses){
  hierstats<-hierarchical.structure.mtDNA.db(ipdb = ipdb_ip,level1 = "sample",level2=h,model="raw",ABGD=F,nperm=1)
  amova_list[[h]]<-hierstats
}
  
#load(file.path(amova_ts_path,"PHIST_TradSpecies_raw_amova_output.Rdata"))
## Summarize AMOVA results
amovastats<-summarize_AMOVA(amova_list,hypotheses,specieslist=unique(ipdb$Genus_species_locus))

WriteXLS(amovastats,ExcelFileName=file.path(amova_ts_path,"PHIST_ts_raw_table_amova_output.Rdata.xlsx"),row.names = T)
save(amova_list,file=file.path(amova_abgd_path,"FST_abgd_raw_amova_output.Rdata"))
save(amovastats,file=file.path(amova_abgd_path,"FST_abgd_table_amova_output.Rdata"))

This takes a long time, so we are not running this code in the document, but loading in the results as we go through each consideration

Visualize AMOVA results

FST and Traditional Species Boundaries

Measure Support

Calculate relative probability from Johnson and Omland 2004 and visualize with a heatmap

minBIC<-sapply(crit_df,function(x){x[which.min(x)]}) 
#J&O box 4 eqn 1. scale() seems to be the way to go here, using the minBIC as the centering vector
crit_df_deltaI<-scale(crit_df, center=minBIC, scale=F) 
#J&O box 4 eqn 4. numerator, plus make it a data frame
crit_df_deltaI_b<-as.data.frame(exp(-0.5*crit_df_deltaI)) 
#J&O box 4 eqn 4. denominator
crit_df_deltaI_sums<-sapply(crit_df_deltaI_b,sum,na.rm=T) 
# this time use the scale argument of scale to divide each column by the corresponding sum
crit_df_relative_prob<-scale(crit_df_deltaI_b,center=F,scale=crit_df_deltaI_sums) 
## Make a heatmap of relative probability of each hypothesis
#melt for ggplot2
relprob<-melt(crit_df_relative_prob)
colnames(relprob)<-c("Hypothesis","Species","Relative_Probability")
#baseplot
rp<-ggplot(relprob,aes(y=Species,x=Hypothesis,fill=Relative_Probability))
#add geom_tile, turn the x-axis elements by 90 degrees, reverse the names on the y-axis, and use a diverging color scheme to highlight hypotheses with >50% rel prob.
rp<-rp+geom_tile()+
  theme(axis.text.x = element_text(angle = 90, hjust = 1))+
  ylim(rev(levels(relprob$Species)))+
  scale_fill_gradient2(low = "blue", mid = "white",
                       high = "red", midpoint = 0.5, space = "rgb",
                       na.value = "grey50", guide = "colourbar")

Non Lab interpolation is deprecated

rp

PHIST and Traditional Species Boundaries

Code is suppressed for the other 3 examples

Measure Support

Calculate relative probability from Johnson and Omland 2004 and visualize with a heatmap

Non Lab interpolation is deprecated

FST and ABGD Boundaries

Measure Support

Calculate relative probability from Johnson and Omland 2004 and visualize with a heatmap

Non Lab interpolation is deprecated

PHIST and ABGD Boundaries

Measure Support

Calculate relative probability from Johnson and Omland 2004 and visualize with a heatmap

Non Lab interpolation is deprecated

Calculate Pairwise Differentiation Stats

By sample

diffstats<-pairwise.structure.mtDNA.db(ipdb=ipdb, gdist = "WC Theta", minseqs = 5, minsamps = 3, mintotalseqs = 0, nrep = 0, num.cores = 2, ABGD = T, regionalization = "sample")


save(diffstats, file="/Users/eric/google_drive/DIPnet_Gait_Lig_Bird/DIPnet_WG4_first_papers/statistics/By_ESU/Pairwise_statistics/DIPnet_structure_sample_WCTheta_ABGD.Rdata")

diffstats<-pairwise.structure.mtDNA.db(ipdb=ipdb, gdist = "Jost D", minseqs = 5, minsamps = 3, mintotalseqs = 0, nrep = 0, num.cores = 2, ABGD = T, regionalization = "sample")


save(diffstats, file="/Users/eric/google_drive/DIPnet_Gait_Lig_Bird/DIPnet_WG4_first_papers/statistics/By_ESU/Pairwise_statistics/DIPnet_structure_sample_JostD_ABGD.Rdata")

Generalized Dissimilarity Modeling

Traditional Species Boundaries

  load("/Users/eric/google_drive/DIPnet_Gait_Lig_Bird/DIPnet_WG4_first_papers/statistics/By_Species/Pairwise_statistics/sample/DIPnet_structure_sample_PhiST.Rdata")
esu_loci <- unique(ipdb$Genus_species_locus)
gdm.full<-list()
solution<-list()
nosolution<-list()
nosolution.full<-list()
nosolution.list<-list()
stats.full<-data.frame(Species_Locus=character(0),WithRegionsDeviance=numeric(0),WithRegionsExplainedDeviance=numeric(0),NoRegionsDeviance=numeric(0),NoRegionsExplainedDeviance=numeric(0),DeltaDeviance=numeric(0),Pvalue_regions=numeric(0),Pvalue_dist=numeric(0),stringsAsFactors = F)
stats<-data.frame(Species_Locus=character(0),Barrier=character(0),WithBarrierDeviance=numeric(0),WithBarrierExplainedDeviance=numeric(0),ImportanceDistanceWithBarrier=numeric(0),ImportanceBarrierWithBarrier=numeric(0),NoBarrierDeviance=numeric(0),NoBarrierExplainedDeviance=numeric(0),ImportanceDistanceWithoutBarrier=numeric(0),DeltaDeviance=numeric(0),Pvalue_barrier=numeric(0),Pvalue_dist=numeric(0),MRDM.rsquared=numeric(0),MRDM.rsquared.pvalue=numeric(0),MRDM.dist.pvalue=numeric(0), MRDM.barrier.pvalue=numeric(0),stringsAsFactors = F)
nostats<-NULL
barriertests<-data.frame(Species=character(0),Region1=character(0),NumPops1=numeric(0),Region2=character(0),Numpops2=numeric(0), Test=logical(0), Solution=logical(0),stringsAsFactors = F)
k<-0
for(gsl in esu_loci){ #gsl<-"Linckia_laevigata_CO1" "Tridacna_crocea_CO1" "Lutjanus_kasmira_CYB" "Acanthaster_planci_CO1"
  
  cat("\n","\n","\n","Now starting", gsl, "\n")
  
  if(any(is.na(diffstats[[gsl]]))){cat("NAs in FST table, No gdm calculated"); next}
  
  if(diffstats[[gsl]]=="Fewer than 3 sampled populations after filtering. No stats calculated"){nostats<-c(nostats,gsl);next}
  
  #pull out the data for this genus-species-locus (gsl)
  sp<-ipdb[which(ipdb$Genus_species_locus==gsl),]
  #clean weird backslashes from names
  sp$locality<-gsub("\"","",sp$locality)
  
  sp$sample<-paste(sp$locality,round(sp$decimalLatitude, digits=0),round(sp$decimalLongitude, digits=0),sep="_")  #sets up a variable that matches the name in Fst table
  sp<-sp[order(sp$sample),]
  # Not all localities are included in Veron's regionalization (e.g. Guam), so get their names and then zap NAs
  nonVeronpops<-unique(sp$sample[is.na(sp$VeronDivis)])
  sp<-sp[!is.na(sp$VeronDivis),]
  
  #subsample Fst 
  gslFST<-diffstats[[gsl]]
  #make a matrix out of gslFST
  gslFSTm<-as.matrix(gslFST)
  
  gslFSTm[which(gslFSTm > 1)] <- 1 #some values that look like 1.0000 are registering as greater than 1
  gslFSTm[which(gslFSTm < 0)] <- 0.0001 #get rid of artifactual negative values
  #gslFSTm<-rescale(gslFSTm)
  #gslFSTm<-gslFSTm/(1-gslFSTm)
  
  #zap weird slashes in the names
  rownames(gslFSTm)<-gsub("\"","",rownames(gslFSTm))
  colnames(gslFSTm)<-rownames(gslFSTm)
  
  #zap the same na populations from the list of non existent pops from VeronDivis
  if(any(rownames(gslFSTm) %in% nonVeronpops)){
    gslFSTm<-gslFSTm[-(which(rownames(gslFSTm) %in% nonVeronpops)),-(which(colnames(gslFSTm) %in% nonVeronpops))]
  }
  
  if(length(rownames(gslFSTm))<5){nostats<-c(nostats,gsl);cat("Fewer than 5 sampled populations");next}
  
  #and filter sp based on the localities that have Fst values
  sp<-sp[sp$sample %in% rownames(gslFSTm),]
  
  #and vice versa
  
  gslFSTm<- gslFSTm[which(rownames(gslFSTm) %in% unique(sp$sample)),which(rownames(gslFSTm) %in% unique(sp$sample))]
  
  #create a locations data frame that has all the localities plus lats and longs and their Veron region.
  locs<-as.data.frame(unique(sp$sample))
  names(locs)<-"sample"
  #locs$Long<-sp$decimalLongitude[which(locs %in% sp$sample)]
  #can't do a unique on sample, lats and longs because some samples have non-unique lats and longs! So I do a join and take the first match.
  locs<-join(locs,sp[c("sample","decimalLongitude","decimalLatitude",
                       "Bowen","Kulbicki_r","REALM")], by="sample", match="first")
  
  
  #sort gslFSTm
  gslFSTm<-gslFSTm[order(rownames(gslFSTm)),order(colnames(gslFSTm))]
  # convert to data frame with popsample names as first column
  gslFSTm<-cbind(sample=locs$sample,as.data.frame(gslFSTm))
  
  ######################################################################
  # 3. Calculate Great Circle Distance
  gcdist_km <- pointDistance(locs[,2:3],lonlat=T)/1000
  #symmetricize the matrix
  gcdist_km[upper.tri(gcdist_km)]<-t(gcdist_km)[upper.tri(gcdist_km)]
  
  #cbind on the sample names
  gcdist_km <- cbind(sample=locs$sample,as.data.frame(gcdist_km))
  
  ####################################################################
  #4. Full Model
  
  cat("Running Full Model")
  
  
  regions<-cbind(locs[,c(1,2,3)],REALM=rep_len(1,length(locs[,1])))
  if(length(unique(locs$REALM)) > 1){
    regions<-cbind(locs[,c(1,2,3)],with(locs, data.frame(model.matrix(~REALM+0))))
  }
  
  gslFSTmll<-cbind(locs[,c(2,3)],gslFSTm)
  
  gdm.format.full<-formatsitepair(bioData=gslFSTm, bioFormat=3, predData=regions,XColumn = "decimalLongitude", YColumn = "decimalLatitude", siteColumn="sample", distPreds=list(gcdist_km))
  
     #zap population pairs with 0 geographical distance between them. Troubling.
    zaps<-which(gdm.format.full$s2.matrix_1==0)
    
    if(length(zaps)>0) {
    gdm.format.full<-gdm.format.full[-zaps,]
    }
  
  # run the full model, and the model with only distance
  fullgdm<-gdm(gdm.format.full)
  distgdm<-gdm(gdm.format.full[,c(1:6,grep("matrix_1",names(gdm.format.full)))])
  
  
  #if no solution, then keep the names of the regions for this species in a matrix with zero values
  if(is.null(fullgdm) | is.null(distgdm)){
        cat("No solution obtained on Full Model")
        regions2<-regions[1,-c(1:3)]
        
        #some futzing to deal with datasets with no predictors but distance matrix
        if(length(regions2)>1){regions2[1,]<-0}
        regions2$matrix_1<-0
        regions2$Species<-gsl
        if(!is.data.frame(regions2)){regions2<-data.frame(matrix_1=0)}
        
        nosolution.full[[gsl]]<-regions2
        nosolution.list[[gsl]]<-gdm.format.full
        plot(gdm.format.full$s2.matrix_1,gdm.format.full$distance, main=gsl)
        next}
  
  # pull out stats
   #difference in deviance
    deltadev.regions<-distgdm$gdmdeviance-fullgdm$gdmdeviance
    
    #percent of null deviance explained by the model with just distance
    explaineddev.dist<-distgdm$explained
    
    gdm.regions.deviance<-fullgdm$gdmdeviance
    gdm.regions.explained<-fullgdm$explained
    gdm.no.regions.deviance<-distgdm$gdmdeviance
    gdm.no.regions.explained<-distgdm$explained
  
  ##############################################################################
    # 4A. Perform Monte-Carlo permutations to develop a null distribution 
    #    of deviance values and determine significance (pvalue)
    gdm.format.rand.full<-gdm.format.full
    rand.deltas.full<-vector() 
    rand.explained.full<-vector()
  
    while(length(rand.deltas.full) < 1000) {
      gdm.format.rand.full$distance<-sample(gdm.format.rand.full$distance,size=length(gdm.format.rand.full$distance))
      gdm.barrier.rand.full<-gdm(gdm.format.rand.full)
      gdm.no.barrier.rand.full<-gdm(gdm.format.rand.full[,c(1:6,grep("matrix_1",
                                                          names(gdm.format.rand.full)))])
    
        # if no solution obtained, go to next gsl
      if(is.null(gdm.barrier.rand.full) | is.null(gdm.no.barrier.rand.full)){next}
      
      deltadev.rand.full<-gdm.no.barrier.rand.full$gdmdeviance-gdm.barrier.rand.full$gdmdeviance
      explained.rand.full<-gdm.no.barrier.rand.full$explained
      
      rand.deltas.full<-c(rand.deltas.full,deltadev.rand.full)
      rand.explained.full<-c(rand.explained.full,explained.rand.full)
    }
    pvalue_regions.full<-length(which(abs(deltadev.regions) < abs(rand.deltas.full)))/length(rand.deltas.full)
    pvalue_dist.full<-length(which(abs(explaineddev.dist) < abs(rand.explained.full)))/length(rand.explained.full)
    
    stats.1<-c(gsl,gdm.regions.deviance,gdm.regions.explained,gdm.no.regions.deviance,gdm.no.regions.explained,deltadev.regions,pvalue_regions.full,pvalue_dist.full)
    
    stats.full[nrow(stats.full)+1,]<-stats.1
    
    gdm.full[[gsl]]<-fullgdm
    
    
}

 
 
 Now starting Acanthaster_planci_CO1 
Running Full Model
 
 
 Now starting Acanthaster_planciNIO_CR 
Running Full Model
 
 
 Now starting Acanthaster_planciPac_CR 
Running Full Model
 
 
 Now starting Acanthaster_planciRedSea_CR 

 
 
 Now starting Acanthaster_planciSIO_CR 
Fewer than 5 sampled populations
 
 
 Now starting Acanthurus_achilles_CYB 
Fewer than 5 sampled populations
 
 
 Now starting Acanthurus_japonicus_CYB 

 
 
 Now starting Acanthurus_leucosternon_CYB 
Fewer than 5 sampled populations
 
 
 Now starting Acanthurus_lineatus_CO1 

 
 
 Now starting Acanthurus_lineatus_CR 
Running Full ModelNo solution obtained on Full Model
 
 
 Now starting Acanthurus_nigricans_A68 

 
 
 Now starting Acanthurus_nigricans_CYB 
Running Full ModelNo solution obtained on Full Model

 
 
 Now starting Acanthurus_nigrofuscus_CO1 
Fewer than 5 sampled populations
 
 
 Now starting Acanthurus_nigrofuscus_CYB 
Running Full Model
 
 
 Now starting Acanthurus_nigroris_CYB 
Running Full ModelNo solution obtained on Full Model

 
 
 Now starting Acanthurus_nigros_CYB 
Fewer than 5 sampled populations
 
 
 Now starting Acanthurus_olivaceus_CYB 
Running Full Model
 
 
 Now starting Acanthurus_sohal_CO1 
NAs in FST table, No gdm calculated
 
 
 Now starting Acanthurus_triostegus_A68 
Running Full Model
 
 
 Now starting Acanthurus_triostegus_CO1 
Fewer than 5 sampled populations
 
 
 Now starting Aipysurus_apraefrontalis_CYB 

 
 
 Now starting Aipysurus_duboisii_CYB 

 
 
 Now starting Aipysurus_eydouxii_CYB 

 
 
 Now starting Aipysurus_foliosquama_CYB 

 
 
 Now starting Aipysurus_fuscus_CYB 

 
 
 Now starting Aipysurus_laevis_CYB 

 
 
 Now starting Aipysurus_mosaicus_CYB 

 
 
 Now starting Aipysurus_pooleorum_CYB 

 
 
 Now starting Amphiprion_bicinctus_CO1 

 
 
 Now starting Amphiprion_omanensis_CO1 

 
 
 Now starting Amphiprion_perideraion_CR 
Running Full Model
 
 
 Now starting Anyperodon_leucogrammicus_CO1 

 
 
 Now starting Apogon_doederleini_CO1 
Fewer than 5 sampled populations
 
 
 Now starting Apogon_doederleini_CR 
Fewer than 5 sampled populations
 
 
 Now starting Arothron_meleagris_A68 

 
 
 Now starting Auxis_rochei_CR 

 
 
 Now starting Auxis_thazard_CR 
Running Full ModelNo solution obtained on Full Model

 
 
 Now starting Caesio_cuning_CR 
Running Full Model
 
 
 Now starting Calotomus_carolinus_A68 

 
 
 Now starting Cantherhines_dumerilii_A68 

 
 
 Now starting Carcharhinus_sorrah_CR 
Running Full Model
 
 
 Now starting Cellana_exarata_CO1 
Running Full ModelNo solution obtained on Full Model

 
 
 Now starting Cellana_sandwicensis_CO1 
Running Full Model
 
 
 Now starting Cellana_talcosa_CO1 
Running Full Model
 
 
 Now starting Centropyge_eibli_CYB 

 
 
 Now starting Centropyge_flavissima_CYB 
Running Full Model
 
 
 Now starting Centropyge_loricula_CYB 
Running Full ModelNo solution obtained on Full Model

 
 
 Now starting Centropyge_vrolikii_CYB 

 
 
 Now starting Cephalopholis_argus_CO1 
Running Full Model
 
 
 Now starting Cephalopholis_argus_CYB 
Running Full Model
 
 
 Now starting Cephalopholis_cyanostigma_CO1 

 
 
 Now starting Cephalopholis_sonnerati_CO1 

 
 
 Now starting Chaetodon_auriga_CO1 
Fewer than 5 sampled populations
 
 
 Now starting Chaetodon_auriga_CYB 
Running Full Model
 
 
 Now starting Chaetodon_collare_CO1 
Fewer than 5 sampled populations
 
 
 Now starting Chaetodon_fremblii_CYB 
Running Full ModelNo solution obtained on Full Model

 
 
 Now starting Chaetodon_leucopleura_CO1 

 
 
 Now starting Chaetodon_lunula_CO1 

 
 
 Now starting Chaetodon_melapterus_CO1 
Running Full ModelNo solution obtained on Full Model

 
 
 Now starting Chaetodon_meyeri_CYB 
Fewer than 5 sampled populations
 
 
 Now starting Chaetodon_miliaris_CYB 
Running Full ModelNo solution obtained on Full Model

 
 
 Now starting Chaetodon_multicinctus_CYB 
Running Full ModelNo solution obtained on Full Model

 
 
 Now starting Chaetodon_ornatissimus_CYB 
Running Full ModelNo solution obtained on Full Model

 
 
 Now starting Chaetodon_trifasciatus_CO1 
Fewer than 5 sampled populations
 
 
 Now starting Chrysoblephus_puniceus_CR 
Fewer than 5 sampled populations
 
 
 Now starting Cirrhitichthys_calliurus_CO1 

 
 
 Now starting Cirrhitichthys_oxycephalus_A68 

 
 
 Now starting Cirrhitichthys_oxycephalus_CO1 

 
 
 Now starting Coryphaena_hippurus_ND1 
Running Full ModelNo solution obtained on Full Model

 
 
 Now starting Cromileptes_altivelis_CO1 

 
 
 Now starting Ctenochaetus_marginatus_A68 

 
 
 Now starting Ctenochaetus_striatus_CR 
Running Full ModelNo solution obtained on Full Model

 
 
 Now starting Ctenochaetus_strigosus_CYB 
Running Full ModelNo solution obtained on Full Model

 
 
 Now starting Culcita_novaeguineae_CO1 
Running Full Model
 
 
 Now starting Culcita_schmideliana_CO1 

 
 
 Now starting Dascyllus_aruanus_CYB 
Running Full Model
 
 
 Now starting Dascyllus_carneus_CO1 

 
 
 Now starting Dascyllus_marginatus_CO1 

 
 
 Now starting Dascyllus_trimaculatus_CO1 
Fewer than 5 sampled populations
 
 
 Now starting Dascyllus_trimaculatus_CR 
Running Full ModelNo solution obtained on Full Model

 
 
 Now starting Dendrochirus_brachypterus_16S 

 
 
 Now starting Dendrochirus_brachypterus_CYB 

 
 
 Now starting Dendrochirus_zebra_16S 

 
 
 Now starting Dendrochirus_zebra_CYB 

 
 
 Now starting Diadema_palmeri_A68 

 
 
 Now starting Diadema_palmeri_CO1 

 
 
 Now starting Diadema_paucispinum_A68 
Fewer than 5 sampled populations
 
 
 Now starting Diadema_paucispinum_CO1 

 
 
 Now starting Diadema_savignyi_A68 
Running Full ModelNo solution obtained on Full Model

 
 
 Now starting Diadema_savignyi_CO1 

 
 
 Now starting Diadema_setosum_A68 
Running Full Model
 
 
 Now starting Diadema_setosum_CO1 

 
 
 Now starting Diodon_holocanthus_A68 

 
 
 Now starting Doryrhamphus_excisus_A68 

 
 
 Now starting Echinothrix_diadema_CYB 

 
 
 Now starting Eleutheronema_tetradactylum_CO1 
Running Full Model
 
 
 Now starting Emydocephalus_annulatus_CYB 

 
 
 Now starting Ephalophis_greyi_CYB 

 
 
 Now starting Epinephelus_areolatus_CO1 
Fewer than 5 sampled populations
 
 
 Now starting Epinephelus_coeruleopunctatus_CO1 

 
 
 Now starting Epinephelus_coioides_CO1 

 
 
 Now starting Epinephelus_fasciatus_CO1 

 
 
 Now starting Epinephelus_flavocaeruleus_CO1 

 
 
 Now starting Epinephelus_longispinis_CO1 

 
 
 Now starting Epinephelus_macrospilos_CO1 

 
 
 Now starting Epinephelus_malabaricus_CO1 

 
 
 Now starting Epinephelus_marginatus_CO1 

 
 
 Now starting Epinephelus_merra_CO1 

 
 
 Now starting Epinephelus_morrhua_CO1 

 
 
 Now starting Epinephelus_multinotatus_CO1 

 
 
 Now starting Epinephelus_ongus_CO1 

 
 
 Now starting Epinephelus_poecilonotus_CO1 

 
 
 Now starting Epinephelus_polyphekadion_CO1 

 
 
 Now starting Epinephelus_spilotoceps_CO1 

 
 
 Now starting Eucidaris_metularia_CO1 

 
 
 Now starting Euthynnus_affinis_CR 
Running Full Model
 
 
 Now starting Forcipiger_flavissimus_A68 
Fewer than 5 sampled populations
 
 
 Now starting Glossogobius_callidus_CYB 
Fewer than 5 sampled populations
 
 
 Now starting Gonodactylellus_viridis_CO1 
Running Full Model
 
 
 Now starting Halichoeres_claudia_CO1 
Fewer than 5 sampled populations
 
 
 Now starting Halichoeres_claudia_CR 
Fewer than 5 sampled populations
 
 
 Now starting Halichoeres_hortulanus_CO1 
Fewer than 5 sampled populations
 
 
 Now starting Halichoeres_hortulanus_CR 
Running Full Model
 
 
 Now starting Halichoeres_hortulanus_CYB 
Fewer than 5 sampled populations
 
 
 Now starting Halichoeres_margaritaceus_CO1 
Fewer than 5 sampled populations
 
 
 Now starting Halichoeres_margaritaceus_CR 
Fewer than 5 sampled populations
 
 
 Now starting Halichoeres_ornatissimus_CO1 
NAs in FST table, No gdm calculated
 
 
 Now starting Halichoeres_ornatissimus_CR 
Running Full ModelNo solution obtained on Full Model

 
 
 Now starting Halichoeres_trimaculatus_CO1 
Running Full ModelNo solution obtained on Full Model

 
 
 Now starting Halichoeres_trimaculatus_CR 
Running Full ModelNo solution obtained on Full Model

 
 
 Now starting Haptosquilla_glyptocercus_CO1 
Running Full Model
 
 
 Now starting Haptosquilla_hamifera_CO1 
Running Full Model
 
 
 Now starting Haptosquilla_pulchella_CO1 
Running Full Model
 
 
 Now starting Heliocidaris_erythrogramma_CO1 

 
 
 Now starting Heliocidaris_tuberculata_CO1 

 
 
 Now starting Heteropriacanthus_cruentatus_A68 

 
 
 Now starting Himantura_leoparda_CO1 

 
 
 Now starting Himantura_tutul_CO1 

 
 
 Now starting Himantura_uarnak_CO1 

 
 
 Now starting Himantura_undulata_CO1 

 
 
 Now starting Holothuria_atra_CO1 
Running Full Model
 
 
 Now starting Holothuria_edulis_CO1 
Running Full Model
 
 
 Now starting Holothuria_whitmaei_CO1 
Running Full ModelNo solution obtained on Full Model

 
 
 Now starting Hoplosquilla_said_CO1 
Running Full Model
 
 
 Now starting Hydrelaps_darwiniensis_CYB 

 
 
 Now starting Hydrophis_atriceps_CYB 

 
 
 Now starting Hydrophis_belcheri_CYB 

 
 
 Now starting Hydrophis_brooki_CYB 

 
 
 Now starting Hydrophis_caerulescens_CYB 

 
 
 Now starting Hydrophis_coggeri_CYB 

 
 
 Now starting Hydrophis_curtus_CYB 

 
 
 Now starting Hydrophis_cyanocinctus_CYB 

 
 
 Now starting Hydrophis_czeblukovi_CYB 

 
 
 Now starting Hydrophis_donaldi_CYB 

 
 
 Now starting Hydrophis_elegans_CYB 

 
 
 Now starting Hydrophis_fasciatus_CYB 

 
 
 Now starting Hydrophis_jerdoni_CYB 

 
 
 Now starting Hydrophis_lamberti_CYB 

 
 
 Now starting Hydrophis_lapemoides_CYB 

 
 
 Now starting Hydrophis_macdowelli_CYB 

 
 
 Now starting Hydrophis_major_CYB 

 
 
 Now starting Hydrophis_obscurus_CYB 

 
 
 Now starting Hydrophis_ocellatus_CYB 

 
 
 Now starting Hydrophis_ornatus_CYB 

 
 
 Now starting Hydrophis_pachycercos_CYB 

 
 
 Now starting Hydrophis_parviceps_CYB 

 
 
 Now starting Hydrophis_peronii_CYB 

 
 
 Now starting Hydrophis_platurus_CYB 

 
 
 Now starting Hydrophis_schistosus_CYB 

 
 
 Now starting Hydrophis_spiralis_CYB 

 
 
 Now starting Hydrophis_stokesii_CYB 

 
 
 Now starting Hydrophis_stricticollis_CYB 

 
 
 Now starting Hydrophis_viperinus_CYB 

 
 
 Now starting Hydrophis_zweifeli_CYB 

 
 
 Now starting Katsuwonus_pelamis_CR 
Running Full ModelNo solution obtained on Full Model

 
 
 Now starting Klyxum_simplex_ND2 

 
 
 Now starting Labroides_bicolor_CR 

 
 
 Now starting Labroides_dimidiatus_CR 
Running Full Model
 
 
 Now starting Labroides_pectoralis_CR 

 
 
 Now starting Linckia_laevigata_CO1 
Running Full Model
 
 
 Now starting Lobophytum_pauciflorum_ND2 

 
 
 Now starting Lutjanus_decussatus_CO1 

 
 
 Now starting Lutjanus_fulvus_CYB 
Running Full Model
 
 
 Now starting Lutjanus_kasmira_CO1 
Fewer than 5 sampled populations
 
 
 Now starting Lutjanus_kasmira_CYB 
Running Full ModelNo solution obtained on Full Model

 
 
 Now starting Microcephalophis_gracilis_CYB 

 
 
 Now starting Mulloidichthys_flavolineatus_A68 

 
 
 Now starting Mulloidichthys_flavolineatus_CR 

 
 
 Now starting Mulloidichthys_mimicus_A68 

 
 
 Now starting Mulloidichthys_mimicus_CR 

 
 
 Now starting Mulloidichthys_pflugeri_A68 

 
 
 Now starting Mulloidichthys_vanicolensis_A68 

 
 
 Now starting Mulloidichthys_vanicolensis_CR 
Running Full ModelNo solution obtained on Full Model

 
 
 Now starting Myripristis_berndti_A68 

 
 
 Now starting Myripristis_earlei_A68 

 
 
 Now starting Naso_brevirostris_CR 
Running Full Model
 
 
 Now starting Naso_caesius_CO1 

 
 
 Now starting Naso_caesius_CR 

 
 
 Now starting Naso_hexacanthus_CO1 
Running Full ModelNo solution obtained on Full Model

 
 
 Now starting Naso_hexacanthus_CR 
Running Full ModelNo solution obtained on Full Model

 
 
 Now starting Naso_unicornis_CR 
Running Full ModelNo solution obtained on Full Model

 
 
 Now starting Negaprion_acutidens_CR 
Fewer than 5 sampled populations
 
 
 Now starting Negaprion_brevirostris_CR 
Fewer than 5 sampled populations
 
 
 Now starting Neoniphon_sammara_CO1 
Fewer than 5 sampled populations
 
 
 Now starting Neoniphon_sammara_CYB 
Fewer than 5 sampled populations
 
 
 Now starting Nerita_albicilla_CO1 
Running Full ModelNo solution obtained on Full Model

 
 
 Now starting Nerita_plicata_CO1 
Running Full ModelNo solution obtained on Full Model

 
 
 Now starting Novaculichthys_taeniourus_A68 

 
 
 Now starting Ostracion_meleagris_A68 

 
 
 Now starting Panulirus_interruptus_CO1 
Fewer than 5 sampled populations
 
 
 Now starting Panulirus_marginatus_CO2 
Running Full ModelNo solution obtained on Full Model

 
 
 Now starting Panulirus_penicillatus_CO1 
Running Full ModelNo solution obtained on Full Model

 
 
 Now starting Parahydrophis_mertoni_CYB 

 
 
 Now starting Parapenaeopsis_coromandelica_CO1 

 
 
 Now starting Penaeus_indicus_CO1 
Running Full ModelNo solution obtained on Full Model

 
 
 Now starting Periclimenes_soror_CO1 
Running Full ModelNo solution obtained on Full Model

 
 
 Now starting Plectorhinchus_chaetodonoides_CO1 

 
 
 Now starting Plectropomus_leopardus_CO1 

 
 
 Now starting Plectropomus_maculatus_CO1 

 
 
 Now starting Polydactylus_macrochir_CR 
Running Full Model
 
 
 Now starting Pomacentrus_coelestis_CO1 
Running Full Model
 
 
 Now starting Pomacentrus_coelestis_CR 
Running Full Model
 
 
 Now starting Pomacentrus_micronesicus_CO1 

 
 
 Now starting Pomacentrus_micronesicus_CR 

 
 
 Now starting Pristipomoides_filamentosus_CYB 
NAs in FST table, No gdm calculated
 
 
 Now starting Protoreaster_nodosus_CO1 
Running Full Model
 
 
 Now starting Pseudoboletia_indiana_CO1 

 
 
 Now starting Pseudoboletia_maculata_CO1 

 
 
 Now starting Pterois_antennata_16S 

 
 
 Now starting Pterois_antennata_CYB 

 
 
 Now starting Pterois_miles_16S 

 
 
 Now starting Pterois_miles_CYB 

 
 
 Now starting Pterois_mombasae_16S 

 
 
 Now starting Pterois_mombasae_CYB 

 
 
 Now starting Pterois_radiata_16S 

 
 
 Now starting Pterois_radiata_CYB 

 
 
 Now starting Pterois_volitans_16S 

 
 
 Now starting Pterois_volitans_CR 
Fewer than 5 sampled populations
 
 
 Now starting Pterois_volitans_CYB 

 
 
 Now starting Pygoplites_diacanthus_CO1 
Fewer than 5 sampled populations
 
 
 Now starting Pygoplites_diacanthus_CYB 
Fewer than 5 sampled populations
 
 
 Now starting Rastrelliger_kanagurta_CR 
Running Full ModelNo solution obtained on Full Model

 
 
 Now starting Rhincodon_typus_CR 
Fewer than 5 sampled populations
 
 
 Now starting Sarcophyton_trocheliophorum_ND2 

 
 
 Now starting Scarus_ghobban_A68 

 
 
 Now starting Scarus_psittacus_CR 
Running Full ModelNo solution obtained on Full Model

 
 
 Now starting Scarus_rubroviolaceus_A68 

 
 
 Now starting Scomberomorus_commerson_CR 
Running Full Model
 
 
 Now starting Scomberomorus_sp_CR 

 
 
 Now starting Sectator_ocyurus_A68 

 
 
 Now starting Sepioteuthis_lessonianaSPA_CO1 
Fewer than 5 sampled populations
 
 
 Now starting Sepioteuthis_lessonianaSPB_CO1 

 
 
 Now starting Sepioteuthis_lessonianaSPC_CO1 
Running Full Model
 
 
 Now starting Siamosquilla_laevicaudata_CO1 
Running Full Model
 
 
 Now starting Siganus_fuscescens_CR 
Running Full Model
 
 
 Now starting Sinularia_densa_ND2 

 
 
 Now starting Sinularia_peculiaris_ND2 

 
 
 Now starting Sphyraena_barracuda_CYB 
Running Full ModelNo solution obtained on Full Model

 
 
 Now starting Sphyrna_lewini_CR 
Running Full Model
 
 
 Now starting Stenella_longirostris_CR 
Running Full Model
 
 
 Now starting Stenella_longirostris_CYB 
Running Full Model
 
 
 Now starting Stenopus_hispidus_CO1 
NAs in FST table, No gdm calculated
 
 
 Now starting Stethojulis_bandanensis_A68 

 
 
 Now starting Thor_amboinensis_CO1 
Running Full ModelNo solution obtained on Full Model

 
 
 Now starting Thunnus_albacares_CO1 

 
 
 Now starting Thyca_crystallina_CO1 
Running Full Model
 
 
 Now starting Tridacna_costata_CO1 
Running Full Model
 
 
 Now starting Tridacna_crocea_16S 

 
 
 Now starting Tridacna_crocea_CO1 
Running Full Model
 
 
 Now starting Tridacna_maxima_16S 

 
 
 Now starting Tridacna_maxima_CO1 
Running Full Model
 
 
 Now starting Tridacna_ningaloo_CO1 

 
 
 Now starting Tridacna_sp_16S 

 
 
 Now starting Tridacna_squamosa_16S 

 
 
 Now starting Tridacna_squamosa_CO1 
Running Full Model
 
 
 Now starting Tripneustes_gratilla_CO1 
NAs in FST table, No gdm calculated
 
 
 Now starting Xenia_sp_ND2 

 
 
 Now starting Zanclus_cornutus_A68 

 
 
 Now starting Zebrasoma_flavescens_CYB 
Running Full ModelNo solution obtained on Full Model

Create a heatmap from GDM regions

dbRDA

# read in the Fst/PhiSt table 
load("~/google_drive/DIPnet_Gait_Lig_Bird/DIPnet_WG4_first_papers/statistics/By_Species/Pairwise_statistics/sample/DIPnet_structure_sample_WCTheta.Rdata")
#load("~/Desktop/DIPnet_structure_sample_PhiST_042817.Rdata")


##1. Dataframe for results
stats<-data.frame(Species_Locus=character(0),constrained.inertia=numeric(0),totalInertia=numeric(0),ProportionConstrained=numeric(0),adj.R2.total=numeric(0),modelF=numeric(0),modelPvalue=numeric(0),pcx_Var=numeric(0),pcx_p=numeric(0),pcy_Var=numeric(0),pcy_p=numeric(0),bestmodel=character(0), constrained.inertia.best=numeric(0), total.inertia.best=numeric(0), proportion.constrained.inertia.best=numeric(0), adj.R2.best.model=numeric(0), stringsAsFactors = F)
stats$Species_Locus<-as.character(stats$Species_Locus)
#stats$Barrier<-as.character(stats$Barrier)
stats$bestmodel<-as.character(stats$bestmodel)

# Make an empty list to save rda output for each species
esu_loci <- unique(ipdb$Genus_species_locus)
all.gsl.rda<-sapply(esu_loci, function(x) NULL)

#make empty data frames to save variance and pvalues
term.vars<-data.frame(gsl=character(0))
term.pvals<-data.frame(gsl=character(0))

###############################################################################
# 2. Subsample for each species of interest, and filter based on Phi_ST table.
for(gsl in esu_loci){ #gsl<-"Linckia_laevigata_CO1" "Tridacna_crocea_CO1" "Lutjanus_kasmira_CYB"
  
  cat("\n","\n","\n","Now starting", gsl, "\n")
  
  if(any(is.na(diffstats[[gsl]]))){cat("NAs in FST table, No dbRDA calculated"); next}
  
  if(diffstats[[gsl]]=="Fewer than 3 sampled populations after filtering. No stats calculated"){all.gsl.rda[[gsl]]<-"Fewer than 5 sampled populations after filtering."; cat("Fewer than 5 sampled populations after filtering.");next}
   sp<-ipdb[which(ipdb$Genus_species_locus==gsl),]
  
  #clean weird backslashes from names
  sp$locality<-gsub("\"","",sp$locality)
  
  sp$sample<-paste(sp$locality,round(sp$decimalLatitude, digits=0),round(sp$decimalLongitude, digits=0),sep="_")  #sets up a variable that matches the name in Fst table
  sp<-sp[order(sp$sample),]
  # Not all localities are included in Veron's regionalization (e.g. Guam), so get their names and then zap NAs
  nonVeronpops<-unique(sp$sample[is.na(sp$VeronDivis)])
  sp<-sp[!is.na(sp$VeronDivis),]
  
  #subsample Fst 
  gslFST<-diffstats[[gsl]]
  #make a matrix out of gslFST, convert negative values to zero
  gslFSTm<-as.matrix(gslFST)
  gslFSTm[which(gslFSTm<0)]<-0.0
  
  #zap weird slashes in the names
  rownames(gslFSTm)<-gsub("\"","",rownames(gslFSTm))
  colnames(gslFSTm)<-rownames(gslFSTm)
  
  #zap the same na populations from the list of non existent pops from VeronDivis
  if(any(rownames(gslFSTm) %in% nonVeronpops)){
    gslFSTm<-gslFSTm[-(which(rownames(gslFSTm) %in% nonVeronpops)),-(which(colnames(gslFSTm) %in% nonVeronpops))]
  }
  if(length(rownames(gslFSTm))<5){all.gsl.rda[[gsl]]<-"Fewer than 5 sampled populations";cat("Fewer than 5 sampled populations");next}
  
  #and filter sp based on the localities that have Fst values
  sp<-sp[sp$sample %in% rownames(gslFSTm),]
  
  #and vice versa
  
  gslFSTm<- gslFSTm[which(rownames(gslFSTm) %in% unique(sp$sample)),which(rownames(gslFSTm) %in% unique(sp$sample))]
  

  if(length(rownames(gslFSTm))<5){all.gsl.rda[[gsl]]<-"Fewer than 5 sampled populations";cat("Fewer than 5 sampled populations");next}
  
  #create a locations data frame that has all the localities plus lats and longs and their Veron region.
  locs<-as.data.frame(unique(sp$sample))
  names(locs)<-"sample"
  #locs$Long<-sp$decimalLongitude[which(locs %in% sp$sample)]
  #can't do a unique on sample, lats and longs because some samples have non-unique lats and longs! So I do a join and take the first match.
  locs<-join(locs,sp[c("sample","decimalLongitude","decimalLatitude"
                       ,"VeronDivis")], by="sample", match="first")
  if (length(unique(locs$VeronDivis)) < 2) {"Only one region!"; cat("Only one region!"); next}
  
  #sort gslFSTm
  gslFSTm<-gslFSTm[order(rownames(gslFSTm)),order(colnames(gslFSTm))]

  
  ######################################################################
  # 4. Calculate Great Circle Distance
  gcdist_km <- pointDistance(locs[,2:3],lonlat=T)/1000
  # and symmetricise it
  gcdist_km[upper.tri(gcdist_km)]<-0
  gcdist_km<-gcdist_km + t(gcdist_km)
   
  ####################################################################### Calculate Overwater Distances#
  #Save for later##

  ##############################################################################
  #5. Create a matrix of regional identities
  
  regions<-with(locs, data.frame(model.matrix(~VeronDivis+0)))   #one of the predictors is superfluous but will get knocked out during RDA
  row.names(regions)<-row.names(locs)
  

    ############################################################################
    # 7. Calculate the principal coordinates
    
    FST.pcoa<-cmdscale(gslFSTm, k=dim(as.matrix(gslFSTm))[1] - 1, eig=TRUE, add=FALSE) #ignore warnings - OK to have negatives according to Anderson
    FST.scores<-FST.pcoa$points
    
    gcdist.pcoa<-cmdscale(gcdist_km, k=2, eig=TRUE, add=FALSE)
    gcdist.scores<-gcdist.pcoa$points
    gcdist.scores<-data.frame("pcx"=gcdist.pcoa$points[,1],"pcy"=gcdist.pcoa$points[,2])
    locs2<-cbind(regions,gcdist.scores)
    
    ###########################################################################
    # 8. Calculate the RDA and extract statistics
    RDA.res<-rda(FST.scores~., data=locs2, scale=TRUE )
    
    #Extract Statistics
    constrained.inertia<-summary(RDA.res)$constr.chi
    total.inertia<-summary(RDA.res)$tot.chi
    proportion.constrained.inertia<-constrained.inertia/total.inertia
    
    adj.R2.total.model<-RsquareAdj(RDA.res)$adj.r.squared
    if(is.na(adj.R2.total.model)){cat("Predictors >= Observations! No solution");all.gsl.rda[[gsl]]<-"Predictors >= Observations! No solution";next}
    
    model.sig<-anova.cca(RDA.res, step=1000)
    modelF<-model.sig$F[1]
    modelPvalue<-model.sig$`Pr(>F)`[1]
    
    terms.sig<-anova(RDA.res, by="term", step=1000)
    pcx_Var<-terms.sig$Variance[1]
    pcx_p<-terms.sig$`Pr(>F)`[1]
    pcy_Var<-terms.sig$Variance[2]
    pcy_p<-terms.sig$`Pr(>F)`[2]
    
    #extract all variance values into a dataframe
    term.var<-as.data.frame(t(terms.sig[,2]))
    names(term.var)<-rownames(terms.sig)
   
    #extract all p-values into a dataframe
    term.pval<-as.data.frame(t(terms.sig[,4]))
    names(term.pval)<-rownames(terms.sig)
    
     
    # scale the variance
    term.var<-term.var/sum(term.var)
    
    # set variance to zero if it is not significant
    term.var[1,which(term.pval[1,] > 0.05 | is.na(term.pval[1,]))]<-0
    
    #set all values to zero if the model itself is not significant
    if(modelPvalue > 0.05){
    term.var[1,which(names(term.var)!="gsl")]<-0
    }
    # add on the species name
    term.pval$gsl<-gsl
    term.var$gsl<-gsl
    
    #merge them on to each dataframe
    term.vars<-merge(term.vars,term.var,all=T)
    term.pvals<-merge(term.pvals,term.pval,all=T)
    
    
    
    #barrier_Var<-terms.sig$Variance[3]   #omitting for now as the number of barriers will differ
    #barrier_p<-terms.sig$`Pr(>F)`[3]
    
    #marg.sig<-anova(RDA.res, by="margin", step=1000)
    #barrier_margVar<-marg.sig$Variance[3]
    #barrier_margp<-marg.sig$`Pr(>F)`[3]
    
    #Model selection
    nullmodel<-rda(FST.scores~1, data=locs2, scale=TRUE )
    forward.model<-ordiR2step(nullmodel, scope=formula(RDA.res), directon="forward", psteps=1000)  #ordiR2step implements Blanchets stopping criterion
    bestmodel<-as.character(forward.model$call[2])
    if (is.null(summary(forward.model)$constr.chi)) {
      constrained.inertia.best<-NA
      total.inertia.best<-summary(forward.model)$tot.chi
      proportion.constrained.inertia.best<-NA
      adj.R2.best.model<-NA
    } else {
      constrained.inertia.best<-summary(forward.model)$constr.chi
      total.inertia.best<-summary(forward.model)$tot.chi
      proportion.constrained.inertia.best<-constrained.inertia.best/total.inertia.best
      adj.R2.best.model<-RsquareAdj(forward.model)$adj.r.squared
    }
     
  
    #save stats
    
    stats_model<-c(gsl,constrained.inertia,total.inertia,proportion.constrained.inertia,adj.R2.total.model,modelF,modelPvalue,pcx_Var,pcx_p,pcy_Var,pcy_p, bestmodel, constrained.inertia.best, total.inertia.best, proportion.constrained.inertia.best, adj.R2.best.model)
    
    stats[nrow(stats)+1,]<-stats_model
    all.gsl.rda[[gsl]]<-RDA.res
    
    
 
  }

write.csv(term.vars,file="./output/dbRDA_term_vars_WCTHETA.csv")
write.csv(term.pvals,file="./output/dbRDA_term_pvals_WCTHETA.csv")

 save(all.gsl.rda,file="./output/multibarrier_dbRDA_WCTHETA.Rdata")
  
write.csv(stats, "./output/multibarrier_dbRDA_JOSTD.csv")
  

Make a heatmap of the variance values